home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / c / newmat.exe / NEWMAT1.CXX < prev    next >
C/C++ Source or Header  |  1991-07-30  |  8KB  |  285 lines
  1. //$$ newmat1.cxx   Utilities and matrix type functions
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5. #include "include.hxx"
  6.  
  7. #include "boolean.hxx"
  8.  
  9. typedef double real;                 // all floating point double
  10.  
  11. #include "newmat.hxx"
  12.  
  13.  
  14. /**************************** error handlers *******************************/
  15.  
  16. void MatrixError(char* message)
  17. {
  18.    cerr << "\nError detected by matrix package: " << message << "\n";
  19.    exit(0);
  20. }
  21.  
  22. void MatrixErrorNoSpace(void* v)
  23. // give error message if v is null
  24. {
  25.    if (v) return;
  26.    cerr << "\nError detected by matrix package: no space\n";
  27.    exit(0);
  28. }
  29.  
  30. /************************* ExeCounter functions *****************************/
  31.  
  32.  
  33. // the next statment may need to be deleted form some compilers
  34.  
  35. int ExeCounter::nreports = 0;
  36.  
  37. ExeCounter::ExeCounter(int xl, int xf) : line(xl), fileid(xf), nexe(0) {}
  38.  
  39. ExeCounter::~ExeCounter()
  40. {
  41.    nreports++;
  42.    cout << nreports << "  " << fileid << "  " << line << "  " << nexe << "\n";
  43. }
  44.  
  45.  
  46. /************************* MatrixType functions *****************************/
  47.  
  48.  
  49. BOOL MatrixType::operator==(const MatrixType& t) const
  50. { return (type == t.type); }
  51.  
  52. BOOL MatrixType::operator!=(const MatrixType& t) const
  53. { return (type != t.type); }
  54.  
  55. MatrixType MatrixType::operator+(const MatrixType& mt) const
  56. {
  57.    Type type1=mt.type;
  58.    if (type==UnSp || type1==UnSp) return UnSp;
  59.    if (type==type1) return type;   // won't work for orthogonal
  60.    if (type==ColV || type1==ColV) return ColV;
  61.    if (type==RowV || type1==RowV) return RowV;
  62.    if (type==EqEl && (type1==Sym || type1==Diag)) return Sym;
  63.    if (type1==EqEl && (type==Sym || type==Diag)) return Sym;
  64.    if (type==EqEl || type1==EqEl) return Rect;
  65.    if (type==Diag) return type1;   // won't work for assym
  66.    if (type1==Diag) return type;
  67.    return Rect;
  68. }
  69.  
  70. MatrixType MatrixType::operator-(const MatrixType& mt) const
  71. {
  72.    return *this+mt;        // won't work for orthogonal or pos def
  73. }
  74.  
  75. MatrixType MatrixType::operator*(const MatrixType& mt) const
  76. {
  77.    Type type1=mt.type;
  78.    if (type==UnSp || type1==UnSp) return UnSp;
  79.    if (type1==ColV) { if (type==RowV) return Diag; else return ColV; }
  80.    if (type==RowV) return RowV;
  81.    if (type==Sym || type1==Sym) return Rect;
  82.    if (type==type1) return type;    // won't work for pos def or assym
  83.    if (type==EqEl || type1==EqEl) return Rect;
  84.    if (type==Diag) return type1;
  85.    if (type1==Diag) return type;
  86.    return Rect;
  87. }
  88.  
  89. BOOL MatrixType::operator>=(const MatrixType& mt) const
  90. {
  91.    Type type1=mt.type;
  92.    if (type==type1) return TRUE;    // by definition
  93.    if (type==Rect || type==RowV || type==ColV) return TRUE;
  94.    if (type==EqEl) return FALSE;
  95.    if (type1==Diag) return TRUE;    // won't work for pos def or assym
  96.    return FALSE;
  97. }
  98.  
  99. MatrixType MatrixType::i() const
  100. {
  101.    switch (type)
  102.    {
  103.    case UnSp:  return UnSp;
  104.    case UT:    return UT;
  105.    case LT:    return LT;
  106.    case Rect:  return Rect;
  107.    case Sym:   return Sym;
  108.    case Diag:  return Diag;
  109.    case RowV:  return Diag;
  110.    case ColV:  return Diag;
  111.    case EqEl:  return Diag;
  112.    case Crout: return UnSp;
  113.    }
  114. }
  115.  
  116. MatrixType MatrixType::t() const
  117. {
  118.    switch (type)
  119.    {
  120.    case UnSp:  return UnSp;
  121.    case UT:    return LT;
  122.    case LT:    return UT;
  123.    case Rect:  return Rect;
  124.    case Sym:   return Sym;
  125.    case Diag:  return Diag;
  126.    case RowV:  return ColV;
  127.    case ColV:  return RowV;
  128.    case EqEl:  return EqEl;
  129.    case Crout: return UnSp;
  130.    }
  131. }
  132.  
  133. MatrixType MatrixType::sub() const
  134. {
  135.    switch (type)
  136.    {
  137.    case UnSp:  return UnSp;
  138.    case UT:    return Rect;
  139.    case LT:    return Rect;
  140.    case Rect:  return Rect;
  141.    case Sym:   return Rect;
  142.    case Diag:  return Rect;
  143.    case RowV:  return RowV;
  144.    case ColV:  return ColV;
  145.    case EqEl:  return EqEl;
  146.    case Crout: return UnSp;
  147.    }
  148. }
  149.  
  150. MatrixType MatrixType::ssub() const
  151. {
  152.    return *this;                  // won't work for selection matrix
  153. }
  154.  
  155. MatrixType::operator char*() const
  156. {
  157.    switch (type)
  158.    {
  159.    case UnSp:  return "UnSp ";
  160.    case UT:    return "UT   ";
  161.    case LT:    return "LT   ";
  162.    case Rect:  return "Rect ";
  163.    case Sym:   return "Sym  ";
  164.    case Diag:  return "Diag ";
  165.    case RowV:  return "RowV ";
  166.    case ColV:  return "ColV ";
  167.    case EqEl:  return "EqEl ";
  168.    case Crout: return "Crout";
  169.    }
  170. }
  171.  
  172. GeneralMatrix* MatrixType::New() const
  173. {
  174.    GeneralMatrix* gm;
  175.    switch (type)
  176.    {
  177.    case UnSp:  MatrixError("Can't build matrix type UnSp");
  178.    case UT:    gm = new UpperTriangularMatrix; break;
  179.    case LT:    gm = new LowerTriangularMatrix; break;
  180.    case Rect:  gm = new Matrix; break;
  181.    case Sym:   gm = new SymmetricMatrix; break;
  182.    case Diag:  gm = new DiagonalMatrix; break;
  183.    case RowV:  gm = new RowVector; break;
  184.    case ColV:  gm = new ColumnVector; break;
  185.    case EqEl:  MatrixError("Can't build matrix type EqEl");
  186.    case Crout: MatrixError("Can't build matrix type Crout");
  187.    }
  188.    MatrixErrorNoSpace(gm);
  189.    gm->Protect(); return gm;
  190. }
  191.  
  192. GeneralMatrix* MatrixType::New(int nr, int nc) const
  193. {
  194.    GeneralMatrix* gm;
  195.    switch (type)
  196.    {
  197.    case UnSp:  MatrixError("Can't build matrix type UnSp");
  198.    case UT:    gm = new UpperTriangularMatrix(nr); break;
  199.    case LT:    gm = new LowerTriangularMatrix(nr); break;
  200.    case Rect:  gm = new Matrix(nr, nc); break;
  201.    case Sym:   gm = new SymmetricMatrix(nr); break;
  202.    case Diag:  gm = new DiagonalMatrix(nr); break;
  203.    case RowV:  if (nr!=1) MatrixError("Illegal Row Vector");
  204.            gm = new RowVector(nc); break;
  205.    case ColV:  if (nc!=1) MatrixError("Illegal Column Vector");
  206.            gm = new ColumnVector(nr); break;
  207.    case EqEl:  MatrixError("Can't build matrix type EqEl");
  208.    case Crout: MatrixError("Can't build matrix type Crout");
  209.    }
  210.    MatrixErrorNoSpace(gm);
  211.    gm->Protect(); return gm;
  212. }
  213.  
  214. void TestTypeAdd()
  215. {
  216.    MatrixType list[] = { MatrixType::UT,
  217.                          MatrixType::LT,
  218.                          MatrixType::Rect,
  219.                          MatrixType::Sym,
  220.                          MatrixType::Diag,
  221.                          MatrixType::RowV,
  222.                          MatrixType::ColV,
  223.              MatrixType::EqEl };
  224.  
  225.    cout << "+     ";
  226.    for (int i=0; i<MatrixType::nTypes(); i++) cout << (char*)list[i] << " ";
  227.    cout << "\n";
  228.    for (i=0; i<MatrixType::nTypes(); i++)
  229.    {
  230.       cout << (char*)(list[i]) << " ";
  231.       for (int j=0; j<MatrixType::nTypes(); j++)
  232.      cout << (char*)(list[j]+list[i]) << " ";
  233.       cout << "\n";
  234.    }
  235.    cout << "\n";
  236. }
  237.  
  238. void TestTypeMult()
  239. {
  240.    MatrixType list[] = { MatrixType::UT,
  241.                          MatrixType::LT,
  242.                          MatrixType::Rect,
  243.                          MatrixType::Sym,
  244.                          MatrixType::Diag,
  245.                          MatrixType::RowV,
  246.                          MatrixType::ColV,
  247.              MatrixType::EqEl };
  248.    cout << "*     ";
  249.    for (int i=0; i<MatrixType::nTypes(); i++)
  250.       cout << (char*)list[i] << " ";
  251.    cout << "\n";
  252.    for (i=0; i<MatrixType::nTypes(); i++)
  253.    {
  254.       cout << (char*)list[i] << " ";
  255.       for (int j=0; j<MatrixType::nTypes(); j++)
  256.      cout << (char*)(list[j]*list[i]) << " ";
  257.       cout << "\n";
  258.    }
  259.    cout << "\n";
  260. }
  261.  
  262. void TestTypeOrder()
  263. {
  264.    MatrixType list[] = { MatrixType::UT,
  265.                          MatrixType::LT,
  266.                          MatrixType::Rect,
  267.                          MatrixType::Sym,
  268.                          MatrixType::Diag,
  269.                          MatrixType::RowV,
  270.                          MatrixType::ColV,
  271.              MatrixType::EqEl };
  272.    cout << ">=    ";
  273.    for (int i = 0; i<MatrixType::nTypes(); i++)
  274.       cout << (char*)list[i] << " ";
  275.    cout << "\n";
  276.    for (i=0; i<MatrixType::nTypes(); i++)
  277.    {
  278.       cout << (char*)list[i] << " ";
  279.       for (int j=0; j<MatrixType::nTypes(); j++)
  280.      cout << ((list[j]>=list[i]) ? "Yes   " : "No    ");
  281.       cout << "\n";
  282.    }
  283.    cout << "\n";
  284. }
  285.